Melting of icosahedral gold nanoclusters from molecular 

dynamics simulations 

Yanting Wang and S. Teitel 
Department of Physics and Astronomy, 
University of Rochester, Rochester, NY 1^627 

Christoph Dellago 
Institute for Experimental Physics, University of Vienna, 
Boltzmanngasse 5, 1090 Vienna, Austria 
(Dated: February 2, 2008) 

Abstract 

Molecular dynamics simulations show that gold clusters with about 600-3000 atoms crystallize 
into a Mackay icosahedron upon cooling from the liquid. A detailed surface analysis shows that 
the facets on the surface of the Mackay icosahedral gold clusters soften but do not premelt below 
the bulk melting temperature. This softening is found to be due to the increasing mobility of 
vertex and edge atoms with temperature, which leads to inter-layer and intra-layer diffusion, and 
a shrinkage of the average facet size, so that the average shape of the cluster is nearly spherical at 
melting. 
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I. INTRODUCTION 



Nanocrystals have quite different physical properties from their corresponding bulk ma- 
terials mainly because of their large surface-to-volume ratio. Among nobel metals, gold 
nanoparticles have already shown their promise for a wide range of applications, such as 
nano-lithographyi, catalysts^, nano-bioelectronic devices^, and ion detection^. Thus knowl- 
edge of the structure and stability of gold nanocrystals is of great importance. 

While bulk gold has an fee crystal structure, the competition between bulk and sur- 
face energies in nanometer sized gold crystallites can result in several different competing 
structures-'-. One such structure which has been observed both in simulationsii^ and in 
experimenta^iiS is the "Mackay icosahedron"^^^, which we will also denote as the "Ih struc- 
ture", consisting of 20 slightly distorted fee tetrahedra, with four {111} facets each, meeting 
at the center to form an icosahedral shaped cluster. The internal facets of the tetrahedra 
meet at strain inducing twin grain boundaries with a local hep structure, leaving the cluster 
with 20 external {111} facets. For an Ih structure with L shells, the magic number of atoms 
needed to construct a perfectly symmetric ideal icosahedron isi^, 

iV^ = y^^ + 5L2 + yL + l. (1) 

In Fig. n] we show the atomic configuration for an ideal Ih structure with the magic number 
of = 2869 atoms (L = 9). Atoms are shaded to indicate local fee, hep, or other structure. 

Different theoretical and numerical models have indicated different limits for the stability 
of such Ih clusters. Using a continuum approach to take into account the strain energy and 
the twin grain boundary energy, Ino predicted that icosahedral clusters should be stable up 
to sizes of 40 000 atomsi^. Similar stability limits for icosahedral gold clusters were predicted 
by Marks using a modified Wulff construction^'^^. More recent atomistic calculation o^^i^^i^^ 
find that, at T = 0, icosahedral gold nanoclusters are the lowest energy structure only 
in a very small size range of tens of atoms; other structures, such as an fee crystal with 
a truncated octahedral shape, become energetically favored at larger sizes. However, in 
molecular dynamics simulations of the freezing of gold nanoparticles Chushak and Bartell^ 
observed icosahedral particles of up to almost 4000 atoms. Simulations by Cleveland et al^ 
found that when fee truncated-octahedral gold clusters with hundreds of atoms were heated, 
they underwent a transformation to the icosahedral structure. 
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FIG. 1: An ideal Ih structure with the magic number of 2869 atoms. Atoms are shaded to indicate 
their local structure: fee is white, hep is gray, and other is black. 

While the theoretical limit of stability of Ih structures thus remains unclear, and their 
formation may involve kinetic rather than strictly equilibrium effecta^'^^'i^'i^*^^^, it is natural 
to suppose that the formation of the Ih structure is related to the very high stability of the 
{111} external surfaces. Simulationsi^ and experiments^^ on bulk slab-like geometries with 
exposed {111} surfaces have shown that, unlike the {100} and {110} surfaces which melt 
below the bulk melting temperature Tm, the {111} surface is a nonmelting surface without 
roughening, wetting, or melting up to the bulk melting temperature T^, and can in fact 
lead to superheating of the solid-". In light of this observation it is interesting to consider 
how the high stability of the {111} facets effects the melting and equilibrium shape of such 
icosahedral clusters. 

In this paper, we will show, by molecular dynamics (MD) simulations, that liquid gold 
clusters with about 600-3000 atoms crystallize into an Ih structure, with a missing central 
atom, upon cooling from the liquid. We then reheat the clusters back through the melting 
transition at temperature T^. Unlike many previous simulations which simulate a continuous 
heating process, we simulate heating in quasi equilibrium, running for a long simulated time 
(43 ns) at each temperature, before increasing the temperature. We pay careful attention 
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to the behavior of the cluster surface, and compute for the first time the average cluster 
shape, as we pass through Tm. Using careful measurements of both bond orientational 
order parameters and atomic diffusion we find that the {111} facets of the cluster surface 
stay ordered and do not premelt or roughen below the cluster melting temperature T^. 
Nevertheless, we find that there is a considerable softening of the cluster surface roughly 
~ 200 K below Tm that can be regarded as due to the "melting" of the atoms on the 
vertices and edges of the cluster. As temperature increases, there is an increasing mobility 
of these atoms leading to intra-layer and inter-layer diffusion, and a shrinkage of the average 
area of the {111} factes. The equilibrium shape progresses from fully faceted, to faceted 
with rounded edges, to nearly spherical just below T^. Throughout this region, the interior 
atoms of the cluster remain essentially perfectly ordered, until is reached. Our results 
refine those of earlier simulations of gold cluster a^i^^i^^ , which used measurements of the 
radial density distribution and the observation of surface diffusion as evidence for a general 
premelting of the cluster surface. A preliminary report of some of our results has been 
presented in Ref.l^- 

II. METHODS 

A. Simulation model and methods 

On modern computers ab initio simulation techniques providing an accurate description of 
interaction energies can be used to simulate systems consisting of up to hundreds of atoma^^. 
However, such methods are too computationally demanding to allow long simulation times. 
Using less accurate but computationally less expensive model potentials, such as the embed- 
ded atom method^^^, the Murrell-Mottram potential^^, the Lennard- Jones potential^, the 
Morse potential^, the many-body Gupta potential^^, or the many-body "glue" potential^, 
one can extend the size of simulated gold nanoclusters to more than ten thousand atoms. 
In this study, we have chosen the many-body "glue" potential because it was found to yield 
an accurate description of the bulk, defect and surface properties of gold^°. In the "glue" 
model, the potential energy of a system of atoms consists of a sum of pair potentials </> 
and a many-body "glue" energy U, 

^ = ^EE</>M + Ef^(^^)- (2) 

i j^i i 
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Here the sums run over all particles, r^j = |rj — rj| is the interatomic distance between 
atoms i and j, and (f){r) is the pair interaction energy. The many-body "glue" energy U{ni) 
depends on an effective coordination number rii of atom i, which is defined by 

j 

Here p(r) is a short-ranged monotonically decreasing function of the interatomic distance r. 
We use the specific forms for the pair potential 0(r), and the glue terms U{n) and p(r), as 



given by Ercolessi et al. in Ref. 
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We will simulate our gold clusters by treating each atom as a classical particle. Newton's 
equations of motion are integrated using the velocity Verlet algorithnt^ with a time step 
of At = 35 fs. Because the many-body "glue" potential uses a cutoff, a cell index method 
can be used to reduce the computational time^^i. In this method, the simulation box is 
divided into cubic cells with side lengths larger than the cutoff distance (3.9 A for the gold 
glue model). When calculating energies and forces one considers only interactions between 
atoms within the same cell and the neighboring 26 cells. This approach reduces the required 
computation time from order A^^ to order A^. On a PC equipped with a 1.5 GHz AMD 
Athlon CPU and 1 GB of memory, 25 000 steps can be carried out per CPU-hour for a 
system of 2624 atoms propagating the system for about 100 ps. 

We will refer to each time step At of the velocity Verlet algorithm as one basic molec- 
ular dynamics (MD) step. Since each MD step is an integration of Newton's equation, it 
necessarily conserves the total energy of the system, and thus simulates a microcanonical 
ensemble. To sample instead in the canonical ensemble, we supplement this basic MD step 
according to two different well known methods. The Andersen thermostali^ is a method 
that mimics a system coupled to a heat bath with constant temperature. At the end of each 
MD step, each particle of the cluster is, with probability p ~ 0.03%, assigned a new velocity 
sampled randomly from the Boltzmann distribution of a given temperature. The Andersen 
thermostat samples both configuration and momentum spaces according to the canonical 
distribution, so that the instantaneous total kinetic energy fluctuates, as is the case for a 
real physical system. However, the Andersen thermostat method does not conserve the total 
linear and angular momenta, and so will cause the system's overall position and orientation 
to drift over the course of the simulation. This would complicate our analysis of cluster 
shape as well as atomic diffusion, since we want to measure these quantities with respect to 



coordinates that stay fixed with respect to the cluster. 

This drawback of the Andersen thermostat can be avoided using the Gaussian isokinetic 
thermostat!^ which keeps the total kinetic energy K of the system fixed at a value cor- 
responding to a given temperature T, K = {3/2)NkBT. We implement this thermostat 
by rescaling all velocities by a constant factor after each basic MD step so as to conserve 
the total kinetic energy. Although the Gaussian isokinetic dynamics does not sample the 
canonical distribution in momentum space, this method does correctly sample configuration 
space and so yields correct equilibrium averages of all structural quantities that depend on 
positional coordinates only. This "constant temperature" method has the advantage of con- 
serving both total linear and angular momenta, thus keeping the cluster at a fixed position 
and orientation. 

In our simulations, our first goal will be to cool our cluster to low temperature from a 
liquid melt to see what solid structure forms. To do this we use the Andersen thermostat 
method since we believe it models more closely the true dynamics of the physical system, 
and hence will incorporate effects that may be due to kinetics rather than just pure ther- 
modynamics. Later, to observe the equilibrium shape and other properties of the cluster, 
we will use the constant temperature MD to heat back through melting. This will keep the 
cluster position and orientation fixed, and so simplify our analysis. 

B. Quantifying structure by bond orientational order parameters 

To determine the crystalline structure of our gold nanoclusters, we will use the method 
of bond orientational order parameter^, which we now review. The idea of the bond 
order parameters is to capture the symmetry of bond orientations regardless of the bond 
lengths. A bond is defined as the vector joining a pair of neighboring atoms. Throughout 
our paper, we will define the neighboring atoms of a given atom i as those atoms which 
have an interatomic distance less than a cutoff radius of 3.8 A, equal to the distance to the 
minimum between the first and the second peaks of the pair correlation function. The local 
order parameters associated with a bond r are the set of numbers. 




(4) 
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where ^(r) and 0(r) are the polar and azimuthal angles of the bond with respect to an 
arbitrary but fixed reference frame, and Yim{d{r),(f){r)) are the usual spherical harmonics. 



Since the bond between atoms i and j may be arbitrarily taken as either rj,- or r 



it is useful to consider only the even-/ bond parameters Qim, since only these are invariant 
to such bond inversions. Global bond order parameters can then be calculated by averaging 
Qimi'r) over all bonds in the cluster, 



1 



Im 



(5) 



^ bonds 

where Nb is the number of bonds. To make the order parameters invariant with respect to 
rotations of the reference frame, the second-order invariants are defined as. 
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and the third-order invariants are defined as. 
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where the coefficients 
quantity. 



mi + 1712 + f^a = 

■ ■) are the Wigner 3j symbols^. It is standard to define a normalized 

Wi 



W, 



3/2 



(8) 



\ Qlm 

which, for a given /, is independent of the magnitudes of the {Qim}- 

The four bond order parameters Q4,, Qq, W4, Wq are generally sufficient to identify 
different crystal structures. In Table|l] we give their values for ideal periodic fee, hep, sc, 
and bcc crystal structures. Since the Ih structure is not periodic, it may in principle have 
values for the bond order parameters that depend on the cluster size. In Fig|2fa) we plot 
the values of the four bond orientational order parameters vs. the cluster size N for several 
ideal Ih structures. Although there is a strong size dependence for small clusters, the bond 
parameters saturate to well defined values as increases, and we list these in Table|l] Note 
that although, for large A^, most of the atoms in the Ih structure have a local fee structure, 
the values of the bond order parameters are different from those of a pure fee crystal; this is 



7 



' 1000 2000 3000 4000 5000 
N 



--Q4 



"'"o 1000 2000 3000 4000 5000 
N 

FIG. 2: The values of (a) the bulk and (b) the surfaee bond orientational order parameters vs. the 
cluster size N, for ideal icosahedral clusters with magic numbers A^. 



TABLE I: Bond order parameters for face-centered-cubic (fee), hexagonal close-packed (hep), sim- 
ple cubic (sc), body-centered-cubic (bcc), liquid, Ih bulk, and Ih surface structures. 



Geometry 


Qa 




W4 


We 


fee 


0.190 94 


0.574 52 


-0.159 32 


-0.01316 


hep 


0.09722 


0.484 76 


0.13410 


-0.012 44 


sc 


0.763 76 


0.353 55 


0.159 32 


0.01316 


bcc 


0.082 02 


0.50083 


0.159 32 


0.01316 


liquid 














Ih bulk 





0.199 61 


-0.159 32 


-0.169 75 


Ih surface 





0.20729 


±0.159 32 


0.169 75 



due to averaging the order parameters over the differing orientations of the 20 fee tetrahedra 
that make up the Ih structure. 

The bond orientational order parameters, averaged over all bonds, will be used to monitor 
global structural changes of our cluster. In particular, a large liquid cluster will have van- 
ishing values for the bond order parameters. Thus the decay of the bond parameters from 
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their finite low temperature values to zero will be a signature of the melting transition. In 
this work we will be particularly concerned with the behavior of atoms on the surface of the 
cluster. Surface structures of nanomaterials can often be quite different from the bulk. We 
therefore will consider separately the bulk bond order parameters, computed by averaging 
over only those bonds connecting atoms that are in the internal "bulk" of the cluster, vs. the 
surface bond order parameters, computed by averaging over only those bonds connecting 
pairs of atoms that lie on the surface of the cluster. In Fig.|2fb) we plot the values of such 
surface bond order parameters for an ideal Ih structure vs. cluster size A^. Again we see 
that they approach well defined constants as increases, except for W4 which oscillates. 
The values of these large A^ surface bond parameters are listed in TableHl Note that for the 
same Ih structure, the bulk bond parameters have different values than the surface bond 
parameters, since they are measuring properties of three dimensional vs. two dimensional 
structures, respectively. 

C. Geometrical analysis of the surface 

A main goal of this work is to quantify the geometrical behavior of the surface of the gold 
nanocluster. In this section we describe the algorithms that we use for this characterization. 

1. Cone algorithm 

The first step in our analysis is to identify the surface particles of the clusters accurately 
according to their geometrical positions. We have developed a new algorithm, which we call 
the "cone" algorithm, to do this. For a given particle, we define an associated cone region 
as the region inside a cone of side length a and angle 6, whose vertex resides on the particle. 
A hollow cone is a cone region with no other particles inside it. We take a particle to be on 
the surface if at least one associated hollow cone can be found. 

The cone algorithm is intrinsically consistent with the general definition of surface parti- 
cles. It can pick up all of the particles on a convex surface. For the particles on a concave 
surface, the precision of identifying a surface atom relies on the choice of the parameters a 
and 6 (see Fig.E)). By visual examination of our generated configurations, we found that 
the two parameters a = 5.0 A and 9 = n/S gave good results for our clusters. For a gold 
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FIG. 3: Schematic of the cone algorithm. 




(a) T=200K 



(b) ] 2U0K 



FIG. 4: Cross section of a gold cluster with 2624 atoms at (a) T = 200 K in an Ih structure and 
(b) T = 1200 K in the liquid, with the 5 topmost layers, as determined by the cone algorithm, 
marked by different gray scales. 



cluster with 5082 atoms, a complete determination of the surface atoms requires less than 2 
s of CPU time. The cone algorithm can also be applied recursively to divide particles into 
surface and sub-surface layers to allow inter-layer analysis. Fig.|^ shows a planar slice cut 
through an instantaneous configuration of a gold cluster with 2624 atoms at (a) T = 200 K 
in an Ih structure and (b) T = 1200 K in the liquid, respectively. From this figure we can 
see that the cone algorithm works well on both solid and liquid configurations. 
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2. Curvature 



Having identified the surface particles, we next want to quantify the surface morphology 
of the clusters. To do this, we compute two different measures of the surface curvature, the 
bond curvature and the maximal local surface curvature, as defined below. 

Our first step is to determine a tangent plane to the cluster surface at each surface particle. 
To do this, we consider the collection of particles determined by the particle of interest and 
all its neighboring particles which are also on the surface. Denote the coordinates of these 
particles by rj = {xi,yi,Zi), i = l...Ns, where A^^ is the number of the particles under 
consideration. We define the tangent plane to pass through the center of mass of these 
particles, and we determine its orientation by minimizing the mean square distance of the 
particles to the plane. Specifically, we solve for this plane as follows. If n = {nx,ny,nz) 
is the unit normal vector of the tangent plane and = {xc, Vc, is the coordinate of the 
center of mass, then the distance from the zth particle to the tangent plane is. 



di 



n 



(9) 



We determine the normal vector n by minimizing X^i df subject to the constraint fi^ 
Introducing an undetermined Lagrange multiplier A, we solve for. 



Eq. (fTUj) leads to the symmetric eigenvalue problem. 



0. 



(10) 



XX -\ XY XZ 
XY YY -\ YZ 
XZ YZ ZZ-\ 



riy 



(11) 



where XY = J^A^i ~Xc){yi — Vc), and other quantities are computed similarly. The smallest 
of the three possible eigenvalues A determines the minimum value of d'^, and its associated 
eigenvector {nx,ny,nz) is the normal vector n of the tangent plane. 

We now define the bond curvature, Cb, of a bond connecting two neighboring surface 
particles. Consider two neighboring surface particles at positions ri and r2 and let fii and 
n2 be the unit normal vectors of their corresponding tangent planes. We can uniquely fit a 
circle that passes through the two points, such that fii and £12 are normal to the circle. The 
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FIG. 5: Schematic for the calculation of the bond curvature. 



bond curvature Cb is then defined in terms of the radius R of this fitted circle, 

1 2sin(^/2) 



Cb = 



R 



(12) 



where 6 = acos(ni ■ £12) is the angle between the two normal vectors and rjj = |ri — is 
the distance between the two particles. The geometry of this fit is illustrated in Fig. El 

Alternatively, we can compute the surface curvature at each particle on the surface as 
follows. Consider a particle on the surface, and all its neighboring particles that are also on 
the surface. Since the geometry of these points roughly defines a surface in 3D space, we 
find the best fit of a paraboloid surface to these points. The maximal local curvature, km, 
and the minimal local curvature, then given by the two principal curvatures of the 

fitted paraboloid. The schematic of this method is shown in Fig.lHl To fit the paraboloid 
at the surface particle Tq, we take the normal vector n to the tangent surface at fq, as 
computed above, and define this to be the local z-axis. We place the origin at tq and then 
choose arbitrary but fixed x and y axes in the tangent plane. We then define the coordinates 
(xj, yi,Zi) of neighboring particle in this coordinate system. We then choose a paraboloid 
function, 

f{x, y) = a^xX^ + 2axyxy + ayyy'^, (13) 
and fit it through the neighboring points by minimizing 

S = Y^{z,-f{x,,y,))\ (14) 

i 

with respect to a^x, o,xy, and ayy. This leads to the following set of linear equations. 



Si -^i Hi X/j -^i Hi X/i -^iUi 



(^xx 

2a. 



\ 



xy 



\ ^yy ) 



( T-x'^z- \ 



(15) 



12 




FIG. 6: Schematic for the calculation of the local surface curvatures. 

Solving Eq. ()15|) for a^x-, cbxyi o-yy, we diagonalize the symmetric matrix with the elements 
a^jy to obtain the two principal axes and the corresponding eigenvalues Ai and A2. The local 
curvatures are then given by ki,2 = 1-^1,2- We will see that the maximal local curvature will 
be very helpful for visualizing the vertex, edge and facet atoms of the cluster surface. 

To test these two methods, we consider the ideal Ih gold cluster with a magic number of 
atoms = 2869, shown in Fig.^ In Fig.[7| (top row) we show the resulting histograms of 
bond curvature and maximal local curvature for this cluster. Both histograms consist solely 
of 5 functions, corresponding to the facets, the edges and the vertices of the Ih structure. For 
the bond curvature histogram there are two separate peaks for the vertices, corresponding 
to the bond curvatures between the vertex atoms and the edge atoms, and the vertex atoms 
and the facet atoms. We also test our method for the average shape of a liquid gold cluster 
with 2624 atoms at T = 1200 K (shown in Fig.|Hl^c)). We show our results in Fig. [7| (bottom 
row). Both histograms have just one peak of finite width, centered at 1/R = 0.047 A~^, 
with i? = 21.5 A the radius of the spherical liquid drop. Note that while the average shape 
of the liquid drop is close to a perfect sphere, the histogram of its curvatures seems relatively 
broad. We have found that this is due to the small discrete number of neighboring surface 
particles that is used to define the fitted paraboloid. Even small deviations of these particles 
from a constant radius can lead to noticeable variations in the fitted curvatures. We find 
that we can reduce the width of the curvature histogram for the liquid cluster by increasing 
the cutoff length used to define neighboring particles, and so have more particles included in 
the fits to the local paraboloids. However, while this improves the calculation for a spherical 
cluster, it makes it worse for high curvature regions near edges and vertices in an Ih cluster, 
where curvature varies rapidly as one moves across the surface. With this understanding of 
our method's limitations, we therefore leave our cutoff as given above. 
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FIG. 7: Histograms of bond curvature and maximal local curvature km of (i) top row: an ideal 
Ih cluster with N = 2869 atoms, and (ii) bottom row: the average shape of a liquid cluster with 
N = 2624 atoms. 

3. Average shape 

At low temperatures, the atoms in the gold cluster remain at well defined equilibrium 
positions and only thermally oscillate around the vicinity of these equilibrium positions. 
The shape of the cluster is thus easily determined from an instantaneous configuration. At 
high temperatures, however, atoms become more mobile and the macroscopic shape of the 
cluster fluctuates dramatically from configuration to configuration. In this case it becomes 
necessary to average over many fiuctuating configurations to define the average cluster shape. 
Since our constant temperature MD conserves total linear and angular momenta and both 
are set to zero, the configurational shape changes we average over represent fiuctuations of 
the surface atoms rather than trivial shifts or rotations of the cluster as a whole. 

At high temperatures, simply averaging the position of each atom throughout the course 
of the simulation docs not give the average shape because the atoms are in general no longer 
bound to specific sites but may diffuse many interatomic spacings through the cluster. We 
therefore use the following approach to define the average cluster shape. We divide the 
surface of the cluster into equal solid angles, and then average the instantaneous surface 
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FIG. 8: A liquid gold cluster with 2624 atoms at T = 1200 K. (a) and (b) are two instantaneous 
configurations, (c) is the shape averaged over 1000 such instantaneous configurations. 

atom positions in each solid angle. This average position in each solid angle then defines 
the profile of the cluster's average shape. This average shape does not contain information 
about the individual surface atom positions, since generally a given solid angle may contain 
the instantaneous positions of different surface atoms at different times. To define our solid 
angle division, we use the best covering spherical codes with icosahedral symmetry^I to 
divide the An total solid angle of the sphere centered at the center of mass into cone cells 
with almost equal solid angles. Choosing different numbers of solid angles results in different 
resolutions; we always choose a number of solid angles that matches as close as possible to 
the number of surface atoms in the cluster. 

We illustrate this method for a liquid gold cluster with 2624 atoms at T = 1200 K. In 
Figs-Efa) and (b) we show two arbitrary instantaneous configurations of the cluster. The- 
oretically, a liquid cluster should have a perfect sphere as the equilibrium shape. However, 
we see that the instantaneous configurations can have noticeably large fiuctuations about 
the average shape. Applying our shape averaging procedure above on 1000 configurations 
sampled at equal time intervals from 43 ns of simulated time, we recover that the average 
shape, shown in Fig.|Hl^c), is a perfect sphere as expected. Note that in Figs.|Hta) and (b), 
the small spheres represent instantaneous atomic positions. In contrast, in Fig|Hl^c), they 
represent not specific atoms, but rather the average surface position within the given solid 
angle. 

D. Atom diffusion analysis 

With enough kinetic energy, atoms can hop around their crystal sites and even travel 
across the whole cluster. The mean squared displacement (MSD)^, Ar^(t), is a convenient 
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way to measure the average movement of a group of atoms. It is defined as 

1 M Ns 

-T^T^Y: [r.fe + t)- v,{t,)f , (16) 

where i — 1 . . . Ns are the positions of the A^^ atoms under consideration, and t is the 
time interval over which the motion takes place. We average over M non-overlapping time 
intervals, with tj = tj-i-\-t. For an infinite three dimensional bulk system, we expect that 
Ar^ — 6Dt as i — > oo, where D is the diffusion coefficient. In a finite cluster, however, the 
MSD will eventually saturate on a length scale comparable to the cluster size. We therefore 
determine the diffusion coefficient D by fitting Ar^(t) to the early time linear part before 
saturation takes place. 

We will also find that a convenient way to visualize individual atomic displacements is 
through an ellipsoid of displacement. We compute this ellipsoid as follows. For a given 
atom traced through K successive configurations for a simulation time t, the mean squared 
displacement correlations are given by the 3x3 matrix C with elements, 

C,.^^j:{r.,-{r,))ir,.-{r.)), (17) 
^ i=i 

where /i,!/ — x, y, z, is the //-coordinate of the atom in configuration i, and (r^u) is the 
average of the coordinate over all K configurations. The probability for the atom to be at 
position r is then approximated as P(r) ~ exp(— |[(r — (r)) • • (r — (r))]), and so the 
surface of our ellipsoid of displacement is given by the equation, 

(r-(r))-C-^-(r-(r)) = l. (18) 

The eigenvectors of C^i, and the square root of their corresponding eigenvalues then define 
the axes and principal radii of the ellipsoid, which we center on the average atom position 
(r) . This ellipsoid provides a convenient visualization of the directional distribution of root 

mean squared displacements over the time t. 

III. RESULTS 

In this section we report on our results. Gold clusters with more than 5000 atoms 
require too much computational time to allow for the long simulation times we want in 
order to explore the equilibrium behavior. Clusters with fewer than a few hundred atoms, 
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however, have large finite size effects due to the larger surface-to- volume ratio. Such smaller 
clusters can undergo transitions between several different crystal structures even at low 
temperature a^i^^i^^ , and they have less sharply defined melting transitions. In this work we 
have therefore simulated several clusters in the range of 600 to 5000 atoms. In our results 
below, we will concentrate on the moderate size of = 2624 atoms, for which we have done 
our most complete and careful analysis. We will also give less detailed results for two smaller 
sizes, = 603 and = 1409, in order to illustrate general trends. Note that these values 
of are not among the magic numbers (see Eq. (^) needed to construct a perfect Mackay 
icosahedron. Nevertheless we will show that these clusters still form Ih structures upon 
cooling. We have also studied several clusters with a magic number of atoms, by explicitly 
constructing the Mackay icosahedron at low temperature, and heating through melting. We 
will give results for sizes^S A^ = 922 and A^ = 5082 in order to compare with the other more 
generic values of A^. 

A. Mackay icosahedra with a missing central atom 

Our initial goal is to cool a liquid cluster through the melting transition to determine 
the ordered structure into which it solidifies. We therefore started with a liquid gold cluster 
with A^ = 2624 atoms which we roughly equilibrated at 1500 K, before cooling to 1200 K 
where we equilibrated longer. We then cooled the cluster down to 200 K, decreasing the 
temperature in intervals of 100 K. At each temperature the system is equilibrated for 5 x 10^ 
steps (21.5 ns) using the Andersen thermostat method. With this cooling method we find 
that our cluster solidifies into an Ih structure^i. 

In Fig. El we show an instantaneous configuration of this A^ = 2624 gold cluster at our 
lowest temperature, T = 200 K. To clarify the geometry of the cluster, we have calculated 
the local curvatures for each surface atom according to the method of Section lTl C 21 and 
in Fig.ini^a) we shade each atom according to the maximal local curvature; the greater the 
curvature, the darker the gray scale. Comparison with Fig.[T] strongly suggests that our 
cluster has an Ih structure. Large curvature regions correspond to edges and vertices, while 
low curvature regions are the fiat {111} facets of the fee tetrahedra. Note that some vertices 
have low curvature; this is because these vertices have their top most atom missing, and so 
form a small locally fiat region. 
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FIG. 9: Ih structure of an = 2624 atom gold cluster at T = 200 K. (a) Surface of an instantaneous 
configuration with atoms shaded according to the maximal local curvature; the larger the curvature, 
the darker the gray scale, (b) The same configuration with the three outer most layers peeled away. 
Atoms are shaded according to their local crystal structure; white is fee, gray is hep, and black is 
"other". 

To further illustrate the Ih nature of our cluster, we have computed the local bond order 
parameters for each atom, averaging over all bonds that connect the given atom to its 
neighbors. Using the values in TableHJ we then identify each atom with its local crystal 
structure. We regard atoms with Q4 > 0.15 and W4 < as having a local fee structure, 
and atoms with Q4 < 0.15 and W4 > as having a local hep structure; all other atoms 
are simply labeled as "other". Because the surface layer and the two sub layers closest to 
the surface exhibit surface reconstruction and have frozen in surface fluctuations, we have 
peeled them away by use of the cone algorithm of Section lTl C 11 The resulting interior of 
the cluster is shown in Fig.El^b), where fee atoms are shaded white, hep atoms gray, and 
"other" atoms black. The Ih structure of the cluster is readily apparent. One clearly sees 
the flat {111} facets of the fee tetrahedra, the edges of the facets corresponding to the hep 
twin grain boundaries, and the vertices with 5-fold symmetry. 

We have also applied the same cooling procedure on smaller gold clusters with = 603 
and N = 1409 atoms. In Fig.^Jwe show the instantaneous configurations of = 603 and 
N = 1409 at T = 200 K, with surface atoms shaded by their maximal local curvature (as 
was done in Fig.Efa) for = 2624). We again clearly see the Ih structure, however for the 
smaller cluster, the edges and facets appear slightly rounded. 

It is interesting to note in Figs.lHlandElthat the fee tetrahedra of our clusters are not 
all of equal size. For a non-magic number of atoms, such as is the case here, this is to 
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FIG. 10: Ih structure of gold clusters with (a) = 603 and (b) N = 1409 atoms at T = 200 
K. The atoms on the surfaces of these instantaneous configurations are shaded according to the 
maximal local curvature; the larger the curvature, the darker the gray scale. 

be expected. However, we have also cooled clusters with magic numbers^ N = 560 and 
= 1414 from liquid to 200 K using the exact same cooling procedure. These clusters also 
formed asymmetric Ih structures with 20 facets of slightly unequal sizes. This suggests that 
our cooling procedure, while slow enough to balance surface vs. bulk free energy and find 
the Ih structures, is not slow enough to achieve the perfect global equilibration which one 
expects would result in perfectly symmetric structures for magic numbers A^. 

An interesting feature of our clusters that can not be seen in Figs.lHlandElis that all of our 
clusters formed with a missing central atom. The energetics of such vacancies at the center 
of Ih clusters were first considered by Boyer and Broughton^° for Lennard- Jones clusters and 
later by Mottet et al.^^ for Cu, Ag, and Au particles. Above a certain material dependent 
critical size the central vacancy lowers the energy of the cluster by partially releasing the 
strain caused by the mismatch of the tetrahedral units. Mottet et al. concluded that for 
gold particles the introduction of the central point defect does not lower the energy enough 
to make the icosahedron competitive with crystallographic octahedra and Wulf polyhedra. 
Their conclusion, however, was based solely on energy calculations which neglect the entropic 
contributions to the free energy at finite temperature. Our finite temperature simulations 
therefore suggest that such a constitutional vacancy can in fact stabilize icosahedral clusters 
of thousands of atoms, making them the observed structure upon cooling. 
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B. Melting and the bond orientational order parameters 

Having determined that our clusters cooled from the liquid have the Ih structure, we then 
heated up the clusters using constant temperature MD instead of the Andersen thermostat, 
so that the total linear and angular momenta are conserved and vanish; this ensures that 
our clusters neither translate nor rotate during the course of our simulations. We heat in 
temperature intervals of 100 K when far from T^, but use smaller intervals when approaching 
Tm. At each temperature the clusters have been equilibrated for 10^ MD steps (4.3 ns), 
followed by 10'' steps (43 ns) to collect data. Our simulation times are more than an order 
of magnitude longer than the ~ 1 ns typically simulated in earlier works^. 

In Fig. ^2 we show the caloric curve (average potential energy per atom vs. temperature) 
for several of our cluster sizes, upon heating. The kink in each curve locates the cluster 
melting transition. Several expected trends^! are clearly seen: (i) the melting temperature 
increases as the cluster size increases, and (ii) the average potential energy per atom increases 
as the cluster size decreases, due to the larger surface-to-volume ratio. No qualitative dif- 
ference is seen between the magic number sizes, N = 922 and 5082, and the others. Note 
that the glue model gives a melting temperature of 1357 K for bulk gold, well above that 
of our biggest cluster—. The experimentally measured melting temperature of bulk gold is 
1337 ¥M. 

We have done our most careful heating for the = 2624 atom cluster, taking fine 
temperature increments near T^. Heating at the above rate of 43 ns per temperature, we 
find that the cluster has a first order melting transition at T = 1075 K. However, when 
we simulated the cluster at the slightly lower temperature of T = 1070 K for more than 
240 ns, we found that it also ultimately melted. Thus the estimates of Tm from Fig.lTTl 
are most likely slightly higher than the true equilibrium values. This superheating that we 
find is perhaps related to the extraordinary stability of the gold {111} surface, as was also 
observed in a slab-like geometry^. 

Next we wish to explore the melting transition from the perspective of the bond orien- 
tational order parameters, defined in Section UTBl We are interested specifically to consider 
the behavior of the surface of the cluster as distinct from the behavior of the interior. We 
therefore use the cone algorithm recursively to group the atoms of the cluster into successive 
layers. The outer most layer of atoms is identified as the surface layer; the atoms immedi- 
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FIG. 11: Caloric curve of Ih gold clusters with = 603, 1409, and 2624 atoms, as well as with 
magic numbers^ of = 922 and 5082 atoms. 



TABLE II: Average numbers of atoms in the surface layer, the sub layers and the bulk of an 
A^ = 2624 atom gold cluster at different temperatures. 



T 


Surface 


Sub layer 1 


Sub layer 2 


Sub layer 3 


Sub layer 4 


Bulk 


400 K 


858.5 ±0.6 


602.8 ±0.8 


428.3 ± 1.1 


307.4 ± 1.1 


207.2 ± 0.8 


219.6 ± 1.1 


600 K 


859.8 ± 1.2 


602.2 ± 1.4 


427.9 ± 1.2 


307.3 ± 1.1 


207.2 ± 0.9 


219.7 ± 1.1 


900 K 


867.7 ±2.4 


594.9 ±2.6 


427.5 ± 1.4 


307.0 ± 1.2 


207.0 ± 1.0 


219.9 ± 1.1 


1060 K 


869.9 ±3.6 


582.4 ±4.0 


436.2 ±3.2 


311.3 ±2.6 


208.6 ± 2.2 


215.7 ±3.2 


1100 K 


874.7 ±3.9 


572.4 ±4.2 


436.2 ± 4.2 


308.7 ±4.0 


209.9 ±3.7 


222.1 ±5.1 



ately below the surface layer are called the first sub layer, then the second sub layer, and 
so on. For the cluster of = 2624 atoms there are a total of 9 such layers. We label the 
atoms lying below the fourth sub layer as "interior" or "bulk" atoms. For A^ = 2624, we 
show in Tableim the number of atoms in each layer for various temperatures up through 
melting. What is immediately apparent is that as the temperature varies within the solid 
phase, T < — 1075 K, the number of atoms in a given layer remains essentially constant, 
within about ~ 5, for all layers below the second sub layer. The surface and top two sub 
layers, however, display a more noticeable variation, suggesting changes on the surface of 
the cluster well below melting. 

Having made this division into layers, we then compute the four bond orientational order 
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parameters Q4, Q^, W4, and Wq, defined in Section|TTH separately for eacli layer and for 
the bulk. In Fig.^J we show our results for the = 2624 atom cluster; Fig. lT^ a) is for 
the interior atoms, while Fig. ll2r b) is for the surface atoms. Comparing to the values listed 
in TablelH or equivalently as shown in Fig.|21 we see that the values we now find at low 
temperature are quite consistent with the bulk and surface values appropriate for an Ih 
structure. The only exception to this is the case of which we find to be approximately 
zero, rather than the negative or positive number shown in Table|l] However we have found 
that, unlike the other bond order parameters, the value of W/^ is extremely sensitive to the 
symmetry of the perfect Ih structure. For deviations from this perfect structure, as is the 
case for our simulated cluster, W4 can vary dramatically. This is evidenced by the very large 
sample to sample fluctuations we found for 1^4, as indicated by the very large error bars 
shown in Fig.[T21 for W4 as compared to the other quantities. We thus conclude that the 
bond orientational order parameters are very consistent with our cluster being a Mackay 
icosahedron. 

In Fig. ll2r a) for the interior atoms, we see that bulk bond orientational order parameters 
remain roughly constant until just above 1000 K, before taking a sharp drop towards zero 
at the same melting temperature, Tm ~ 1075 K, as found from the caloric curve of Fig.lTTl 
Thus the bond orientational order parameters give a good signature of the melting transition. 
The sharp decrease of the bond parameters indicates that the interior atoms remain with 
a highly ordered Ih structure until just before melting. Note that the values in the liquid 
above are not identically zero, but have small finite values due to the finite size of the 
liquid cluster; this effect is biggest for Q^. 

In Fig. ll2r b) for the surface atoms, we again see that the bond orientational order pa- 
rameters remain with their Ih values at low temperatures, and then vanish towards zero at 
the same as for the bulk atoms. Thus we reach one of our most important conclusions: 
the presence of finite surface bond orientational order up until the bulk melting transition 
indicates that the surface {111} facets of the Ih structure do not premelt, but rather the 
surface facets melt at the same temperature as the bulk. The absence of any sharp features in 
the surface bond orientational order parameters below Tm suggests that there are no other 
types of surface phase transitions below T^. There is, however, one noticeable difference in 
the behavior of the surface bond orientational order parameters as compared to the bulk. 
We see that the surface parameter Qq starts a noticeable decrease from its low temperature 
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FIG. 12: Bond orientational order parameters of the N = 2624 atom cluster for (a) the interior 
atoms, and (b) the surface atoms. Sample error bars, representing configuration to configuration 
fluctuations, are shown. 

value at T ~ 800 K, considerably below the melting T^. We interpret this as a softening of 
the surface, and we will present the reason for this behavior in the following section. 

We have also measured the bond orientational order parameters for the first through 
fourth sub layers of the cluster. Since Qq is the bond parameter that most clearly shows the 
surface softening, in Fig.Hniwe plot the value of Qe vs. temperature, for surface, interior, 
and each of the four sub layers. Since each layer has a slightly different value of Qe at low 
temperature, we plot the normalized values Qq{T)/Qq{AOOK) so as to better compare their 
relative behaviors. We see that both the surface and the first sub layer show almost identical 
softening as Tm is approached. However all the deeper sub layers show almost identical 
behavior as the interior atoms, with almost no softening until Tm. Thus the softening 
phenomenon is seen to be largely confined to the top two layers of the cluster and does not 
propagate more deeply as is approached; below the top two layers, the cluster remains 
almost as ordered as at low temperatures, until just before melting. 

We have also tested the sensitivity of our definition of the "interior" atoms of the cluster, 
by redefining it to be all the atoms below the surface layer. However, as might be expected 
from Fig.El computing the bulk bond orientational order parameters defined this way gives 
no qualitative change from the behavior seen in Fig. lT^ a). 

In Figs. El to El we show similar plots of interior and surface bond orientational order 
parameters for our other cluster sizes, N = 603, 922, 1409 and 5082. We see the same 
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FIG. 13: Normalized bond orientational order parameter Qq{T) /Qq{400K) for the surface, interior, 
and various sub layers of the = 2624 atom cluster. 
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FIG. 14: Bond orientational order parameters of the N = 603 atom cluster for (a) the interior 
atoms, and (b) the surface atoms. 

qualitative behaviors as in Fig.^J with surface and bulk melting at the same temperature. 
This melting temperature, which increases with cluster size, agrees with the values found 
from the caloric curves of Fig.^2 Surface softening tracks the melting temperature and 
starts to be noticeable about 200 K below Tm. The surface softening is somewhat enhanced 
for the smaller cluster sizes. There appears to be no qualitative differences for our magic 
numbei^ clusters, = 922 and 5082, as compared to the other sizes. 
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FIG. 15: Bond orientational order parameters of the magic number = 922 atom cluster for (a) 
the interior atoms, and (b) the surface atoms. 
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FIG. 16: Bond orientational order parameters of the = 1409 atom cluster for (a) the interior 
atoms, and (b) the surface atoms. 

C. Average shape and surface curvature 



To understand the physical manifestations of the surface softening that is indicated by 
the surface bond orientational order parameters, we now look at the average shapes of our 
cluster, computed according to the method of Section lTl C 31 We focus first on our cluster 
of = 2624 atoms. For this case, we have divided the Air total solid angle into 842 



37. We have chosen this 



almost equal solid angles, using the icosahedral covering of Ref. 
number since it corresponds as close as possible to the typical number of surface atoms 
in the cluster (see TableUT)). At each temperature over 1000 instantaneous configurations. 
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FIG. 17: Bond orientational order parameters of the magic number = 5082 atom cluster for (a) 
the interior atoms, and (b) the surface atoms. 

sampled at equal intervals throughout the simulated time of 43 ns, have been included in 
our average. We show the resulting average shapes in Fig.^1 We present results for the 
following temperatures: 400 K, representing the low temperature configuration in which 
thermal fluctuations are negligible; 600 K where one starts to notice small changes in the 
surface; 900 K where substantial softening of the surface orientational order parameter Qq 
is observed; 1060 K, just below the melting ^ 1075 K; and T = 1100 K just above T^. 

In the top row of Fig.^1 we show pictures constructed similarly to that in Fig.|Hfc) of 
Section lTl C 3\ which showed the average shape of the hquid cluster at T = 1200 K. The 
small spheres represent the average position of the surface within the given solid angle. 
Additionally, we have now shaded these spheres according to the value of the maximal local 
surface curvature, as we did earlier in Fig.Efa) for an instantaneous configuration at T = 200 
K; the darker the gray scale, the larger the curvature. This method of shading is used to 
highlight any edges and facets that are on the cluster surface. The view point for these 
pictures is taken at infinity, so as to show a full hemisphere of solid angle. 

In the bottom row of Fig.^J we show the corresponding average shapes using a smooth 
3D contour plot with overhead lighting. The view point for these bottom row pictures is 
now taken to be a finite distance from the cluster, in order to highlight the straight edges 
and 5-fold symmetry about the vertices. 

The pictures in Fig. ^1 illustrate the following scenario as the cluster is heated. At low 
temperatures the cluster is almost fully faceted, with flat facets meeting at sharp edges and 
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FIG. 18: Average shapes of an iV = 2624 atom cluster at 400, 600, 900, 1060 and 1100 K. The 
top row shows each of the discretized sohd angles of the surface shaded according to the value of 
the maximal local curvature; the darker the gray scale, the larger the curvature. The viewpoint of 
these pictures is set to infinity, to show a full hemisphere of solid angle. The bottom row is the 
corresponding smooth contour plot, with a finite viewpoint so as to highlight the straight edges 
and 5-fold symmetry about the vertices. 

vertices. By 900 K the facets have shrunken slightly in size and the edges and vertices have 
noticeably rounded. At 1060 K, just below melting, the facets have shrunken to almost 
negligible size, and the cluster is almost spherical. Above melting, the cluster is essentially 
a perfect sphere. 

As a way to quantify the cluster shapes we have computed the bond curvatures Cb and 
the maximal local surface curvatures km, as defined in Section lll C 21 In Figs. lTW aVfd) we 
show histograms of bond curvature Cb for the four temperatures 600, 900, 1060 and 1100 
K. The solid curves show the histograms of bond curvatures as computed over the surface 
bonds of the average cluster shape shown in Fig.^1 In contrast, the dashed curves show 
the histograms of bond curvatures computed for an instantaneous cluster configuration, and 
then averaged over the 1000 instantaneous configurations saved in our simulated time of 43 
ns. Note that for the histograms for the average shape, where since we are dealing with only 
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one average configuration we liave relatively few points in our histogram, we have smoothed 
our data using a Gaussian smoothing function with a width of 4 bins. The bin size here 



to a region where the surface is locally concave; an example of when this can happen is near 
a vertex which is missing its top most atom. 

Both Figs. El and 1201 illustrate the same scenario. Consider first the histograms of the 
average cluster shapes. At low temperatures, the histograms show a strong peak at zero, 
representing the low curvatures of the large flat facets. The histograms also show either a 
second peak or plateau at higher curvature, with a long high curvature tail, representing 
the large curvatures at edges and vertices. We can compare these results against those in 
Fig. [71 for the ideal Ih structure. In the liquid above Tm ^ 1075 K, the histograms have a 
single sharp peak at finite curvature, representing the uniform curvature of the spherical 
liquid cluster. Just below melting, at T = 1060 K, the histograms similarly show a single 
peak near that of the liquid, only noticeably broader than for the liquid; this indicates the 
shrinkage of the flat facets to negligible size and a rounded cluster that is not yet a perfect 
sphere. 

Comparing the histograms for the average vs. the instantaneous shapes, the latter are 
in general broader, most especially for the liquid cluster. This demonstrates the presence of 
strong thermal shape fluctuations about the average shape. Particularly interesting are the 
histograms for 1060 K, just below Tm, in Figs. imT c) andEHTc). and for 1100 K, just above Tm, 
in Figs. fmT d) andEOfd). The histograms for the average shape are symmetric Gaussian like 
peaks about an average curvature, corresponding to a spherical or nearly spherical cluster. 
The histograms for the instantaneous configurations, however, are skewed in shape. They 
have a low curvature peak and a broad high curvature tail, somewhat similar to what is seen 
at lower temperatures. This suggests that the instantaneous configurations can still develop 
small local facets on the surface. A similar observation has previously been made by Lewis et 
al. for smaller clusters^Si. Fluctuations of the edges and vertices of these local facets lead to 
an effective diffusion of the facet upon the cluster surface; averaging over these fluctuations 
results in a smoothing out of the facets to negligible size when one considers the average, 
rather than the instantaneous, cluster shape. We have seen evidence for this scenario by 
visual inspection of instantaneous cluster configurations. For example, in the instantaneous 
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FIG. 19: Histograms of bond curvature Cb for the average cluster shape (sohd hues) and the 
instantaneous cluster configurations (dashed lines) at (a) T = 600K, (b) T = 900 K, (c) T = 1060 
K, and (d) T = 1100 K. The cluster size is = 2624 atoms. 

liquid cluster configuration shown in Fig.|Hfb) one clearly sees a flat edge along the bottom. 

For comparison with other sizes, we show in Fig.|^ average cluster shapes for our = 
1409 atom cluster at temperatures 800 K and 900 K, where — 925 K. In Fig.|221we show 
average shapes for our magic number— N = 5082 atom cluster at temperatures 1000 K and 
1140 K, where now Tm — 1150 K. The gray scale in these figures is the same as used in 
Fig.^1 Again we see facets shrinking, and the cluster becoming more spherical, as is 
approached. 



D. Diffusion of atoms 



In this section we present further evidence that the physical mechanism behind the surface 
softening is indeed the diffusion of atoms on the vertices and edges of the cluster. We will 
consider in this section only the cluster of = 2624 atoms. 

We start by first considering the mter-layer mixing of atoms in the cluster, defining an 
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FIG. 20: Histograms of maximal surface curvature km of the average cluster shape (solid lines) 
and the instantaneous cluster configurations (dashed lines) at (a) T = 600K, (b) T = 900 K, (c) 
T = 1060 K, and (d) T = 1100 K. The cluster size is TV = 2624 atoms. 




T=80QK 



FIG. 21: Average cluster shapes for an N = 1409 atom cluster at temperatures 800 K and 900 K, 
where ~ 925 K. 

inter-layer mixing parameter (n) as follows. At each temperature we label the atoms in the 
initial configuration by an integer n' = 0, 1, 2, . . . , 5 according to whether the atom is on the 
surface, in the first sub layer, second sub layer,. . . , or interior of the cluster. At the end of the 
simulation for that temperature, we assign a new integer n to each atom, according to which 
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FIG. 22: Average cluster shapes for an = 5082 atom cluster at temperatures 1000 K and 1140 
K, where Tm ^ 1150 K. 
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FIG. 23: Interlayer mixing parameter (n) vs. T, for atoms initially on the surface, in the first sub 
layer, . . . , and in the interior. The cluster size is N = 2624 atoms. 

layer the atom is now in. In Fig.|2niwe plot (n) vs. T, where n is averaged separately over 
each group of atoms indexed by their initial layer number n'. When (n) differs noticeably 
from the initial n', it indicates significant inter-layer mixing of the atoms from layer n' into 
other layers. From Fig.|2Slwe see that noticeable inter-layer mixing takes place between the 
surface and the first sub layer as low as 700 K; these two layers are almost evenly mixed by 
950 K, more than 100 K below the melting Tm ~ 1075 K. As Tm is approached, additional 
layers start to mix together. At Tm and above, all layers are evenly mixed during the course 
of the simulation, indicating that in the liquid all atoms diffuse equally throughout the entire 
cluster. 

Next we consider the diffusion of the atoms in the cluster by computing the mean squared 
displacements, Ar^(t), defined in Eq. (fT?)|) . We compute Ar^(t) separately for each layer of 
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FIG. 24: Mean squared displacements for the N = 2624 atom cluster averaged over the atoms in 
the surface layer, first through fourth sub layers, and interior for (a) 600 K, (b) 900 K, (c) 1060 K, 
and (d) 1100 K. 

the cluster (and the interior) by averaging only over the atoms that are initially in a given 
layer. In Figs. l24r a)-(d) we plot our results for Ar'^{t) vs. t, layer by layer, for the four 
different temperatures, 600, 900, 1060 and 1100 K. Note that since atoms in different layers 
can mix (see Fig.l23p. the division into different layers in Fig.|^ contains some ambiguity; an 
atom initially in the first sub layer, for example, might during the course of the simulation 
wind up on the surface, however we continue to average its motion in with that of the first 
sub layer. 

Several expected features are apparent in Fig.|m In Fig.|2H^d) at 1100 K, above the 
melting — 1075 K, we see that all layers behave roughly the same, saturating at Ar^ ~ 
600 A^. This is consistent with a liquid cluster of radius ~ 21 A, in which all the atoms can 
diffuse throughout the entire cluster, no matter which layer they were initially in. At the 
low temperature of 600 K, where the average cluster shape remains almost fully faceted, the 
results in Fig.|2Ha) show that diffusion is almost negligible. Even for the two top layers. 
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FIG. 25: Diffusion coefficients D vs. T for different layers of the N = 2624 atom cluster. The inset 
shows an expanded range for D in a temperature range below melting, 700 - 1050 K. 

atoms on average move less than one inter-atomic spacing (~ 3 A) over the observation 
time of 20 ns. At 900 K, where the edges and vertices of the average cluster shape have 
noticeably rounded, we see in Fig. E^ b) that diffusion in the top two layers is significant, 
with atoms on average traveling a root mean square distance equal to several inter-atomic 
spacings. The second sub layer also shows a noticeable diffusion but all more inward atoms 
diffuse negligibly. At 1060 K, just below the melting Tm ~ 1075 K, we see that all atoms 
are diffusing a significant amount throughout the cluster, with the top two layers almost 
reaching the long time saturation value ~ 600 found in the liquid. 

In Fig.|2Hlwe plot the diffusion constant D vs. T for each of the cluster layers, obtained 
by fitting to the early time linear part of the curves in Fig.|211 If we fit our diffusion 
constant for the surface layer to the simple form D = Doexpl—Ejs^/k^T) to extract the 
activation energy, = —d{lnD)/d{l/kBT), we find the values of = 0.21 eV at low 
temperatures, ~ 500 K, where the cluster is fully faceted. At high temperature, ~ 1200 
K, in the liquid, we find Ep^ = 0.35 eV. Note that the first value corresponds to surface 
diffusion, while the second value corresponds to bulk diffusion in the liquid (sine once the 
cluster has melted, atoms initially on the surface easily diffuse into the bulk). To compare 
with previous simulations, Boisvert et al.— did a first principles calculation for the gold 
{111} surface at low temperatures and found Ea = 0.22 ± 0.03 eV, in good agreement with 
our value. Chushak and Bartell^ reported the value of E^ = 0.25 eV using the EAM model 
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for a liquid gold nanocluster. Considering the tendency of the EAM model to systematically 
give lower energy values (as pointed out in Ref.l44). our result for the liquid is in reasonable 
agreement. 

A seeming paradox concerning our diffusion results of Fig.|21]is that at low temperatures 
the diffusion of atoms in the first sub layer appears to be greater than that for atoms on the 
surface. This can be explained by noting that the Ar^ in Fig. [^represent an average over 
all the atoms in a given layer. As we will show below, at low temperatures, atoms along the 
edges and on the vertices of a given layer are more mobile than a typical atom in that layer. 
Since the fraction of such edge and vertex atoms is larger in the first sub surface layer than 
on the surface, atoms in this layer have a larger average mobility. When the temperature 
increases to 1060 K, most of the atoms in the two layers are now diffusing, and the average 
mean squared displacements Ar^ of the two layers become roughly equal. 

A more serious issue is how to reconcile our results of Fig.|^ showing noticeable surface 
diffusion below Tm, with our claim that the surface {111} facets remain ordered and do 
not premelt below Tm, as indicated from the finite values of the bond orientational order 
parameters shown in Fig. lT^ b). One possibility is that the surface layer does in fact melt 
at a well defined temperature below T^, but that orientational order in the liquid surface 
is maintained due to the presence of an effective periodic substrate formed by the ordered 
sub layers below the surface. However we do not believe that this is the case. Even if 
orientational order in a liquid surface were preserved by the presence of the ordered sub 
layers, one would still expect to see some kink or other feature in the bond orientational 
order parameters at the surface melting transition. In contrast, we find in Fig. 112b that the 
bond order parameters go smoothly, though near steeply, to zero. Instead of the above 
scenario, we believe that the surface diffusion that we observe below Tm is due not to atoms 
on the {111} facets, but rather due to the atoms along the vertices and edges of the surface. 
As temperature increases, the facets shrink in size, the edges get rounder and broader, and 
the effective number of such diffusing atoms increases. Just below Tm the facets have shrunk 
to almost negligible size, the bond order parameters have decreased to a corresponding small 
but finite value, and most of the atoms on the surface are now diffusing. 

To estimate the number of atoms in each layer that are diffusing, we use the following 
criterion. We compute the number of atoms in each layer that have moved a distance 
of more than 8 A within 20 ns of simulated time. The cutoff of 8 A is chosen since it 
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in 20 ns, for the = 2624 atom cluster. 

is the distance between the third and fourth peak of the pair correlation function, thus 
representing a distance roughly between the third and fourth nearest neighbor; we assume 
that an atom which can move this far is in fact diffusing, rather than just undergoing thermal 
motion about a fixed average position. We find our results to be qualitatively insensitive to 
choosing a smaller cutoff length of 6 A (the distance between the second and third peaks 
of the pair correlation function). If Fig. |^ we plot the fraction of such "moved" atoms vs. 
temperature, for the surface, sub layers, and interior of the cluster. We see that only the 
surface and the first sub layer have a significant fraction of moved atoms below melting. 
This fraction steadily increases with temperature, and approaches unity at T ~ 1000 K, just 
below Tin. We interpret the unmoved fraction as those atoms on the ordered {111} facets, 
which shrink in size as T approaches Tm. Close enough to T^, when the facets become so 
small that they are only a few atoms across, it becomes easy for atoms on or near the edge 
of a facet to exchange with mobile atoms in the surrounding "liquid" of edge atoms; hence 
even such facet atoms can ultimately diffuse throughout the cluster, and the fraction of 
moved atoms can approach unity below T^. Indeed the concept of liquid vs. solid become 
somewhat ambiguous when referring to such small surface the {111} facets just 

below Tm. 

Combining all the above we infer the following scenario for diffusion at low temperatures: 
only atoms on the surface and in the first sub layer show any noticeable diffusion well below 
Tm. The atoms in these two layers that diffuse are the same atoms which mix between 
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the two layers (see Fig. 1221), these are the atoms along the edges and vertices of each 
layer. The atoms from the first sub layer diffuse by migrating first to the surface, and then 
diffusing upon the surface, until mixing back into the first sub layer. As the temperature 
increases to T^, the facets shrink and the number of diffusing edge atoms increases, until all 
surface atoms are diffusing just below T^. 

To substantiate the above picture, we plot in Fig.|2Z|the displacement ellipsoids, defined 
in Section lllDl for all atoms initially on the surface of our = 2624 cluster. We show results 
for temperatures 400, 600, 900, 1060 and 1100 K, corresponding to the same temperatures 
for which we showed the average cluster shape in Fig.^1 In the top row we show ellipsoids 
averaged over a simulated time of 1.075 ns. We expect that atoms which are diffusing, with 
Ar^ ~ t, should have their displacement ellipsoid roughly double in size when the time 
interval goes up by a factor of four. Hence in the bottom row of Fig. EH we then show results 
for a simulated time of 4.3 ns, i.e. four times longer than the top row. We observe the 
following. At 400 K there is no observable diffusion of surface atoms. At 600 K we see 
diffusion of atoms at the vertices of the Ih cluster. At 900 K we see stronger diffusion at 
the vertices, as well as diffusion along the edges. One also can see several of the ellipsoids 
oriented normal to the surface, indicating atoms which are mixing in with the first sub layer. 
Atoms at the centers of the facets remain without diffusion. At 1060 K and above most of 
the atoms are clearly diffusing. 



IV. DISCUSSION AND CONCLUSIONS 

We have carried out long time equilibrium molecular dynamics simulations to study the 
behavior of gold nanoclusters cooled from the liquid, and their subsequent melting upon 
reheating. For three different generic cluster sizes, A^ = 603, 1409, and 2624, we found that 
the cooled clusters formed a slightly asymmetric Mackay icosahedral (Ih) structure with a 
missing central atom. 

Using the above clusters cooled from the melt, as well as several other "magic number" 
Mackay icosahedra with up to A^ = 5082 atoms that we constructed by hand^^ at low 
temperature, we slowly heated these clusters up through melting. Measuring surface and 
bulk bond orientational order parameters we find a sharp cluster melting transition at a 
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t = WE T = 600K T = 900K T = 1060K T=1100K 




FIG. 27: Ellipsoids of displacement at 400, 600, 900, 1060 and 1100 K for the cluster of TV = 
2624 atoms. Each ellipsoid is centered at the average position of the given atom, and shows the 
directional distribution of root mean squared displacements. The top row gives results obtained 
for a simulated time of 1.075 ns, while the bottom row is for 4.3 ns. 

temperature Tm(iV) that increases with cluster size, and that the {111} facets on the surface 
do not prcmclt but remain ordered up until Tm. The surface bond parameters however 
decrease from their perfect Ih values significantly below T, indicating a softening of the 
surface prior to melting. We find that the onset of this surface softening appears to track 
the size dependent melting, occuring roughly 200 K below Tia{N). 

Looking at the average shape of our clusters, we see that this surface softening corresponds 
to a rounding of the edges and vertices of the cluster, with a corresponding shrinkage of the 
{111} facet area. Just below the melting Tm the average cluster shape is nearly spherical. 
As the temperature increases towards melting, and in the liquid above T^, instantaneous 
cluster configurations can display large thermal fiuctuations about this average shape. 

Measuring the diffusion of atoms in the cluster, we conclude that the mechanism for this 
surface softening is the onset of diffusion of atoms at first the vertices and then the edges 
of the cluster surface, as temperature is increased. As temperature further increases, the 
mobility of these atoms increases, and more and more atoms near the edges of the facets 
participate in this diffusion, until the number of atoms remaining stationary on the facets 
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becomes almost negligibly small near T^. Simultaneous with this increasing diffusion is an 
increase in interlayer mixing, with surface and first sub layers mixing first, and then deeper 
layers mixing in as one approaches close to T^. 

A similar rounding of edges and shrinking of facets occur in the theory of the equilibrium 
shape of macroscopically large crystals, where the continuum Wulff construction^^ can be 
applied. In this theory, the shrinkage of facets is associated with approaching a roughening 
transition of the faceted surface, and the facet length shrinks proportionally to the inverse 
of the roughening correlation length^. We do not believe that this theory explains the 
results for our clusters. Firstly, it is generally believe d^^i^^ that the {111} gold surface that 
forms the facets of our cluster does not have a roughening transition below the bulk melting 
transition. This is consistent with our observation that the surface softening of our clusters 
seems to track the size dependent cluster melting temperature rather than approaching a 
size independent onset temperature as would be expected if there was a true thermodynamic 
roughening transition. Moreover, the vanishing of facets at the roughening transition occurs 
within the context of the crystalline state; no diffusion of atoms need be involved. In our 
case it is clear that diffusion of atoms along the vertices and edges plays an important role. 
We therefore believe that the phenomena we observe in our simulated nanoclusters are due 
specifically to the finite, relatively small, size of our clusters, for which continuum approaches 
are not valid and one must consider the atomistic nature of the system. The rounding of 
edges and shrinkage of facets that we observe are better attributed to a "melting" of the 
cluster edges, which then spreads out into the ordered facets as the temperature increases 
towards melting. 
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